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Abstract. In the most popular approach to the numerical study of the Anderson metal- 
insulator transition the transfer matrix method is combined with finite-size scaling ideas. 
This approach requires large computer resources to overcome the statistical fluctuations and to 
accumulate data for a sufficient range of different values of disorder or energy. In this paper 
we present an alternative approach in which the basic transfer matrix is extended to calculate 
the derivative with respect to disorder. By so doing we are able to concentrate on a single value 
of energy or disorder and, potentially, to calculate the critical behaviour much more efficiently 
and independently of the assumed range of the critical regime. We present some initial results 
which illustrate both the advantages and the drawbacks of the method. 
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1. Introduction 

The Anderson (1958) transition in three-dimensional electronic systems has been extensively 
studied over recent years (Kramer & MacKinnon 1993). It is well known that increasing the 
fluctuations of the random potential in three-dimensional systems (3D) causes a transition 
from metal to insulator. This transition occurs due to a change in the nature of the electron 
wavefunctions under the influence of the disorder. In the insulating phase the magnitude 
of the potential fluctuations is sufficiently large to localise the wavefunctions to a region 
of space. As the disorder is reduced below a critical value the electron states become 
delocalised allowing conduction through the system. In models of disordered systems such as 
the Anderson Hamiltonian this change can be measured numerically by observing the change 
in the decay length of the transmission probability (the correlation length £) using the transfer 
matrix method (MacKinnon & Kramer 1981, MacKinnon 1994). Since the development 
(Abrahams et al. 1979) of a scaling theory for the zero-temperature dc conductance of 
disordered electronic systems much numerical work has focused on calculating the critical 
exponent for the correlation length (MacKinnon & Kramer 1981, MacKinnon 1994, Slevin & 
Ohtsuki 1997, Slevin & Ohtsuki 1999). 

The assumption of one-parameter scaling implies that the renormalised length Ai = 
£i/M should follow the scaling equation: 

dlnAi 

dlnM = * l(lnAl) (1) 
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with solutions of the form A 1 = fiiM/^). When plotted against disorder, W, for different 
sizes, M, the curves, fx, all cross at a fixed point Xi = corresponding to the metal-insulator 
transition. £ x is the correlation length for the infinite system which diverges close to the 
critical point, W c , as ^ ~| W — W c \~ u . Near the critical point the scaling equation can be 
linearised: 

Ai =A hc + a(W-W c )M^. (2) 

The standard method of calculating the critical exponent v is to analyse the finite size scaling 
in terms of this equation. However, taking the logarithm of the derivative with respect to 
disorder gives us: 




In a + - In M. (3) 

v 



This expression provides a way of calculating the critical exponent from the derivative of the 
correlation length of the finite system Ai. In previous work (MacKinnon 1990) it has been 
noted that the calculation of the derivative dAi/dVK would provide useful complimentary 
information about the scaling and a more direct approach to calculating the critical exponent. 
However direct estimation of the derivative has up to now been plagued by large numerical 
errors associated with numerical division when using a finite difference approximation 
[Ai(W + AW) - Ai(W)]/AW. 

Our new method is an adaptation of the standard transfer matrix approach to calculate the 
derivative, T, of the transfer matrix, T, from which it is possible to calculate the expectation 
value of dAi/dW, and hence find v. This approach is applicable for 3D cubic systems 
(M=L=6, 8, 10,12,14) rather than for quasi-lD systems. In contrast to quasi-lD systems, 
where the conductance is described by only one correlation length, the conductance of cubic 
systems is given by a sum over contribution from a series of correlation lengths. Here we 
present preliminary results for the scaling of each correlation length individually as well as 
for the conductance as a whole. 



2. Computational Method 

2.1. The traditional transfer matrix 

The Anderson (1958) model for non-interacting electrons on a simple cubic lattice in zero 
magnetic field is described by the Hamiltonian: 

H = J2ti\i)(i\+Y, V \i)(3l (4) 

i i=£j 

where \i) is the atomic orbital at site i, and V = 1 for nearest neighbours and zero otherwise. 
The disorder is introduced by allowing the on-site energies q to behave stochastically. The 
site energies are chosen according to the probability distribution 

p{e) = W- 1 6{\W -\e\), (5) 

where 9 is the step function so that p(e) is a box distribution. The Schrodinger equation 
H(j) = E(f) for a system with cross section M x M and length L can be rewritten in recursive 
form: 
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where the 2M 2 x 2M 2 matrix, T, is the so-called the transfer matrix and the M 2 x M 2 
matrix, Hj is the Hamiltonian of a single slice disconnected from the rest of the lattice. The 
application of Tj to a slice gives the solution of the Schrodinger equation for the next slice. In 
this way the complete wavefunction for a disordered sample of length L can be generated by 
repeated application of the transfer matrix: 



^0 



(7) 



oo) 



The standard transfer matrix method then consists of calculating the asymptotic (i.e. L 
behaviour of the eigenvalues of the Hermitian matrix Q: 

Q = T f T (8) 

by repeated orthogonalisation of the columns of T (MacKinnon & Kramer 1983). 
Alternatively, for a finite L, a system in Landauer (1970) geometry, i.e. embedded between 
perfectly ordered leads, may be represented in the form 

t 



Q = U f 



II T « 



AITITA 



u 



(9) 



where the columns (rows) of U (U T A) are the right (left) eigenvectors which diagonalise T in 
the absence of disorder and 

1 
1 



(10) 



The correlation lengths from which the scaling properties are calculated, are related to the 
eigenvectors \zi > of Q by 

Q|zj >= exp(zj)|zj > where z t = 2ctiL = (11) 

for i = 1,2, ...2M 2 . In the quasi ID limit the on are called Lyapunov exponents. 

The normalised eigenvectors for the system with zero disorder can be grouped together 
to form a 2Af 2 x 2M 2 matrix, U, in which first M 2 columns are made from the 2M 2 x M 2 
sub-matrix (U + ) of eigenvectors representing waves travelling in the positive z-direction, 
and the last M 2 columns from those travelling in the opposite direction, (U_). The set of 
left-handed eigenvectors is given by U T A. The eigenvectors can be used to transform the 
real-space transfer matrix (6) into the wave representation: 

T w = U T AT L U, (12) 

where is the real-space transfer matrix for a system of length L. The transfer matrix 
in the wave representation, T w , can be expressed in terms of the transmission and reflection 
matrices: 



t - r'(t / ) _1 r r'(t')- 1 
-(t'r 1 ' (IT 1 
This suggests that r' and t' may be calculated from the matrix products: 



i-i 



t 



r't' 



= ( 1 ) T„ 



(1 ) T„ 



U T AT L U -U T A[]T n U 

n=l 
L 

u t + a n t « u 



U|AT L U 



(13) 



(14a) 



(14ft) 



n=l 
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2.2. Stabilisation of the Iteration 

In practice the repeated multiplications of the matrix U by successive transfer matrices gives 
rise to numerical instabilities that require attention. The matrix gradually becomes dominated 
by the largest eigenvalues and the smaller ones are lost owing to the finite accuracy of the 
numerical process. It is these eigenvalues that are needed when the inverse is taken to calculate 
t'. To avoid losing them the algorithm must be modified. After about ten iterations, and well 
before the small eigenvalues are lost, the relevant information is 'stored' by multiplying from 
the right by a stabilising matrix, Y: 
10 

[]T„U Y = ZY (15) 

n=l 

The matrix Y may be chosen to be the inverse of the top half of the Z (or the upper triangular 
matrix which orthonormalises the columns of Z): 

zy =(V) y =(gy)- < 16 > 

This process of multiplying by the inverse matrix is repeated approximately every ten 
iterations. The product of the Y is stored separately and is used at the end to solve for the 
reflection/transmission coefficients. Applying this method to equations (14a) and (14b) results 



in the following relations: 

t' _1 Y = A, (17a) 

rt' X Y = B, (17/7) 
where the matrices A, B and Y are known. The numerically stable solutions are: 

t' = YA \ (18a) 

r =BA \ (lSb) 



2.3. Calculating the derivatives 



Our new approach extends the transfer matrix method to calculate T = dT/dVF from which 
the d£, / dW can then be found. To preserve the multiplicative behaviour a larger matrix K is 
used: 







n 



t, 







where Tj has the simple form: 
t, = 



H, 


where Hi is the derivative of the Hamiltonian of slice i: 

/e'l ... 

^ ... 



H, 



dH 

dW 




V o 



£ M 2 -1 










£ M 2 



(19) 



(20) 



(21) 
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where is a random number between —0.5 and 0.5 from (5). The larger size of the matrix 
increases the computational time required but because of the additional information contained 
in the derivative only one calculation per system size is needed to calculate v, rather than many 
calculations over a wide range of disorders as in the standard method. The computational 
details of how to calculate the derivatives of g and & from Tj are described in section 2. In 
section 3.1 preliminary results of the scaling of these quantities is presented. 

This method cannot be applied to long systems (i.e. quasi-lD) because, from Osledec's 
theorem, we would obtain a limiting Kl matrix of eigenvalues which cannot be related back 
to the eigenvalues of T. A possible solution would be to stabilise T and T simultaneously 
so that they would both lead to limiting matrices of eigenvalues. However, two different 
sets of vectors are required to orthogonalise T and T and this creates problems at the end of 
the calculation when it becomes necessary to remove the stored stabilisation vectors. The 
only way to overcome these difficulties is to consider small cubic systems which do not need 
stabilising. In this work cubic systems of size M = 6, 8, 10, 12 were used. 

In the polar decomposition the transfer matrix can be parameterised as (Mello et al. 
1988): 

T = UTV, (22) 




■t 



(23) 



where Ui are arbitrary M x M unitary matrices and A is a real diagonal matrix with positive 
elements If time-reversal symmetry holds then U3 = uf and U4 = U^. In this 

parameterisation the transmission matrix t' can be written as: 

t' = U 3 (l + A)- 1 / 2 U 2 (24) 

and the conductance as: 

g = Tr(t't't) = Tr(u 3 (l + A) _1 uJ), (25) 

M 1 

(26) 



= E 



, 1 + A< 

The Lyapunov exponents for the ensemble are derived by substitution of the polar 
decomposition (23) in the matrix product Q (8,9) 



/ Ui( l + 2A)uI 2u lv /A(l + A)u 3 \ 
V \2uy\(l + X)u\ l4(l + 2A)u 3 J { ' 

Q is Hermitian positive so the eigenvalues are positive real numbers and the flux conservation 
constraint, Q 1 = AQA requires that the eigenvalues come in inverse pairs. Strictly speaking 
the Lyapunov exponents on (11) are defined only in the quasi- ID limit, but in this work 
we will also use the term for arbitrary L. Applying Osledec's theorem of random matrix 
products (Pichard & Sarma 1981a, Pichard & Sarma 1981/?), the Lyapunov exponents are 
self-averaging for large L. Hence, taking the limit of large L is statistically equivalent to 
an ensemble average over many different realisations of the disorder. Computationally, for 
samples of width M close to the transition, convergence of the first Lyapunov exponent to 1% 
accuracy is achieved for a length L ~ 10 4 x M. The conductance can be expressed in terms 
of the Lyapunov exponents (Pichard 1984): 

N 1 

■9 = 2 E 1 + cosh(2Laj) < 28 > 

i=l 
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In the quasi-lD limit (L 3> 1) the conductance can be written as: 

N / N \ 

g ~ e ~ 2a * L = e ~ 2aiL 1 + S e~ 2( "*- Ql) (29) 

i=l \ i=2 / 

where the A" ~ M d_1 Lyapunov exponents are ordered and ol\ is the smallest. As the length 
of the system is increased the fractional contribution to the current from the higher Lyapunov 
exponents rapidly decreases. The conductance is dominated by the first Lyapunov exponent 
and the system is similar to the ID case where there is only one Lyapunov exponent (one 
conductance channel) and all the states are localised. Since there are no other relevant length 
scales present the localisation length is equivalent to the correlation length of the finite system, 
£i, and is defined as: 

i_ = - lim ^B. = lim ai (30) 

From the Hellman-Feynman theorem the expectation value of the derivative of the 
eigenvalues, n, of t'^t' can be calculated: 

dn t dt*f 

dW =u >MV Ui (31fl) 



AW 

where Ui is the eigenvector corresponding to n. To find the Ui and r, it is necessary to 
perform an initial calculation of the transfer matrix as outlined in section 2.1. Then Ui is used 
in a calculation of the derivative from equation (31b). The expectation value of the derivative 
is only calculated for the smallest ten Lyapunov exponents to avoid the necessity of stabilising 
the matrices. If the expectation value of the whole matrix was to be calculated, the smallest 
expectation values would rapidly become insignificant compared to the largest, and would be 
lost in the rounding errors. In general the expectation value for a quasi ID system cannot be 
calculated because as L increases the result rapidly becomes independent of initial conditions, 
i.e. independent of the initial vector fa. 

Starting from equation (31b) we derive the form used in the computational calculations. 
From (14a) we have: 

t' 1 = U T AT L U (32) 

where U_ correspond to eigenvectors in the leads with negative k z . Taking derivatives and 
using equation (19) we obtain: 

from the identity t' 1 t' = 1 we can write: 

^ = -t'^t' (34) 
AW dW { V 

therefore we have: 

= -t'U T At L U_t' (35) 

AW 

where the transmission matrices come from the initial calculation and still contain the 
evanescent modes. Only at this point can we remove the evanescent modes (indicated by 
application of Q) and calculate the expectation value: 

^(1 + A,) -1 = -0 4 Q W A ( ? )n iC *( o ) U - t,Q(/)4 (36) 
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where Q is an M x M' matrix with all elements zero, except for Q nn = 1, 1 < n < N 
which correspond to the N non-evanescent modes. The actual computation starts with the 
initial vector on the right and by successive matrix products moves through to the left. The 
derivative of A 4 can be calculated (Pichard 1984): 

TTTF = - 4 ^[ cosh_1 ( 2A ' + !)) 2 sinh(cosh- 1 (2A i + l))]" 1 (37) 
dW dW 

These values were found to be in agreement with the derivative as calculated from the normal 
method, w (A 17 5 — A 16 5 ). The critical exponent is calculated by fitting a straight line to 
equation (3). 

The derivative of the conductance g can be calculated from the derivative of the LE's 
(37). The one-parameter scaling theory is formulated in terms of the conductance and the 
scaling behaviour of this quantity is considered first in the next section before considering 
each LE individually. 



3. Cubic systems 




0.28 1 ■ 1 1 1 ■ 1 ■ 1 ■ 1 1 1 

16.4 16.6 16.8 17.0 17.2 17.4 17.6 

W 

Figure 1. Conductance, g in units of e 2 /h, vs. disorder, W, for M = 6 ( ), 8 ( ), 

10 (---), 12 ( )&14( ). 



3.1. Scaling of the conductance 

To estimate the critical point we have plotted (figure 1) the conductance g against disorder 
W for system sizes M = 6, 8, 10, 12, 14. The value of critical disorder, W c , is not expected 
to depend on the quantity being considered. Therefore, a size-independent critical point at 
approximately W — 16.5 would be anticipated to be in agreement with the values calculated 
from both energy statistics (Hofstetter & Schreiber 1993, Zharekeshev & Kramer 1994) of 
cubes and from the analysis of the 1st LE in quasi-lD systems. However, it is clear that the 
apparent critical disorder is above W = 16.5 within the error bars of the calculation. The 
results are also qualitatively the same for log(<?). This movement of the critical point from 
the expected value of W — 16.5 is probably due to the presence of the leads connected to the 
ends of the cube. This is in agreement with other results (Slevin et al. 2000, Kawarabayashi 
et al. 2000) which find that the conductance is size dependent at W = 16.54 with periodic 
boundary conditions. However, with hard wall boundary conditions the conductance is found 
to be virtually size-independent at W — 16.54 (Slevin & Ohtsuki 1997). In energy statistics 
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calculations the boundary conditions affect the critical distributions but do not appear to shift 
the critical point (Braun & Pascaud 1998). 

Table 1. At disorder W = 16.5 the number of iterations, open modes and the 1st LE for cubic 
samples , E = 0. 



M 


iterations 


modes 


Ai 


6 


100000 


21 


0.671 ±0.0015 


8 


50000 


45 


0.694 ±0.002 


10 


30000 


61 


0.709 ±0.003 


12 


20000 


95 


0.718 ±0.003 


14 


5000 


123 


0.729 ±0.007 



-0.4 - 



ln(D) 




-2.4 - 



-2.8 1 — 1 — 1 — 1 — ' — 1 — 1 — 1 — ' — 1 — 1 — 1 — ' — 1 — 1 — 1 — ' — 1 — 1 — 1 — 1 
1.7 1.9 2.1 2.3 2.5 2.7 

ln(M) 

Figure 2. In D vs. In M for W = 16.5 and M = 6,8, 10, 12, 14. The gradient is equal to 

1/v. D = (dlog( 5 )/dW0 ( ), (dlog( S )/dW0 / <log( S )> ( ), (dg/dW) / (g) 

(---), (dg/dW) ( ). 

Applying equation (3) with g as the scaling parameter gives us: 

ln(D) = b+ ilnM (38) 

with D — dg/dW\ w . If log g is used as the scaling parameter then D = 1/g dg/dW\ w . 
Also considered were D — (dg/dW\ w )/{g) and D = (d\og{g) / dW\ w )/(\og(g)). The 
plot of equation (38) is given in figure 2 and the critical exponent found by a 'least-squares- 
fit' of the gradient is given in table 2. The error in the value of g is considerably smaller 
than the error in the gradient and can be ignored when calculating the standard deviation 
of the fit. From the error bars on the graph we can see that log(g) produces slightly more 
accurate results because it decreases the influence of the infrequent but very large values of 
the derivative. The M = 14 data is not used in the calculation of the critical exponents. 
The results from the same calculation performed at W c — 17.5 are also given in table 2. On 
average the critical exponents appear to be slightly lower in value for W c — 17.5. Clearly the 
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Table 2. Critical exponents calculated from the conductance of cubes with sizes M 
6,8,10,12 

W D v x 2 Q 

16.5 1.85 ±0.12 0.7 0.7 

dW 

16.5 ( d '^ g) )/( lo g(3)) 1-34 ±0.066 1.0 0.6 

16.5 1.24 ±0.095 0.3 0.9 
dW 

16.5 ^.^L^/(g) 1.65 ±0.17 0.2 0.9 

17.5 dl °g(3) 1.34 ±0.073 1.4 0.5 
dW 

17-5 ( dl °^ 9) )/<log(g)) 1.30 ±0.068 1.2 0.5 

17.5 — 1.23 ±0.14 0.4 0.8 



17.5 ^^/( fl ) 1.28 ±0.15 0.3 0.9 



observed shift of the critical point is completely independent of our new method to calculate 
the critical exponent. However this deviation from scaling theory may well be prejudicing 
the critical exponent obtained from the new method. Nevertheless there is general agreement 
with previous calculations of the critical exponent which give a value of v m 1.5. 

There are thus two deviations from the expected results, both of which are probably 
due to finite size effects: the shift in the apparent critical disorder and the dependence of the 
calculated exponent on the particular average of g used. However, another possible factor is 
the properties of the leads for which the number of open channels represents an additional 
degree of freedom. This requires further investigation. 



3.2. Scaling of the Lyapunov exponents 

For cubic systems attached to leads the contribution of the 1st 'Lyapunov exponent' § to 
the value of the conductance in the critical region can be calculated. At W = 16.5 
the 1st LE makes up 85% of the conductance and approximately 74% of the derivative 
of the conductance. Therefore the higher Lyapunov exponents do have physical meaning 
even though the dominant contribution is from the 1st exponent. Scaling higher Lyapunov 
exponents is analogous to scaling energy statistics at large energy scales (such as P(i, s) for 
large i) (Taylor 1998), since both correspond to short length scale events. 

The 1st LE is plotted against disorder in figure 3. Within the error bars of the calculation 
it can be concluded that the critical point is above W — 16.5 as was the case with the 
conductance g. This shift of the critical point is even clearer in the more accurate data for 
the 6th LE (figure 3). Qualitatively similar results occur for (Aj). 

In figure 4 equation (3) is plotted with D = (dA~ 1 /dW) for W = 16.5, i.e. the first 

§ The Lyapunov exponents on are strictly defined only in the quasi-lD limit but in this work the definition is 
extended to include arbitrary L (11). 




Figure 3. Dependence of (Ai) (upper figure) and (Ag) (lower figure) on disorder W for cubic 

systems of size M. Data points at W = 16.5, 17.0, 17.5. M = 6 ( ), 8 ( ), 10 

(---), 12 ( ),&14( ). 
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-1.2 - 
ln(D) _! 4 . 

-1.6 - 




-1.8 1 1 1 1 1 1 1 1 1 1 ' 

1 .7 1.9 2.1 2.3 2.5 2.7 

ln(M) 



ln(D) 




_-\ 2 ' 1 1 1 1 1 1 1 1 1 1 

'1.7 1.9 2.1 2.3 2.5 2.7 

ln(M) 

Figure 4. Scaling of the Lyapunov exponents. Slope is equal to 1/u. D = (dAT 1 /dW^ at 
i = l (upper figure), i = 6 (lower figure). 

LE, (A -1 ), is being scaled. 

From one parameter scaling the gradient of this line is equal to the inverse of the critical 
exponent. 

The accuracy of the data increases when we consider higher LE's and the result for the 
6th LE is shown in figure 4. In theory this should be a straight line but the accuracy of the 
data is high enough to observe a definite curvature indicating deviations from one-parameter 
scaling theory. Although the curve can be fitted by a straight line the quality-of-fit measure 
X 2 is too large (and Q -C .1) which indicates the curvature cannot be ignored. At W = 17.5, 
which is nearer the apparent critical point (figure 5), the results are qualitatively the same and 
the curvature in the \n(D) vs ln(M) plot is again observed for the higher Lyapunov exponents. 
This type of curvature indicates that scaling corrections exist for the gradient. A full analysis 
in terms of this equation was not possible because the number of data points is too small to 
obtain sufficient accuracy. In future studies with a full compliment of data from M = 4 to 
14 it is possible that the functional form of the corrections could be found and subtracted to 
obtain a crossing at W — 16.5. It is impossible to tell if the curvature also exists for the 1st 
LE because it could be masked by the larger error bars of the data. 

The results of the straight line fit to calculate the critical exponent are given in figure 6. 
As mentioned, although the critical exponents calculated from the higher Lyapunov exponents 
are more accurate the x 2 becomes too large. However the fits illustrate that the scaling 
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Figure 6. Critical exponents calculated from 1/Aj (higher curve) and from Aj (lower curve) 
for i = 1 — 6. W = 16.5 (upper figure), W = 17.5 (lower figure) 
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behaviour is approximately in agreement with previously calculated values of the critical 
exponent. 

4. Conclusions 

We have introduced a new approach to calculating the critical behaviour of the Anderson 
transition in disordered systems, which has the potential to overcome some of the numerical 
difficulties which have plagued such calculations in the past. As the method is only requires 
a calculation at a single value of disorder it eliminates the difficulties associated with the 
necessity of fitting to a range of values of disorder and the associated uncertainties associated 
with the determination of the range of disorder over which the critical behaviour is valid. The 
initial results are promising, although the necessity of working with cubes rather than long 
systems introduces deviations from the expected behaviour. We expect these deviations to 
be explicable in terms of the usual corrections to scaling (Slevin & Ohtsuki 1997, Slevin & 
Ohtsuki 1999, Ohtsuki et al. 1999, Slevin & Ohtsuki 2001). 
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